% This script estimates peer, contextual and individual effects
% Date: 2019-03-14
clear; close all; clc;
% filepath = 'C:\Users\xt9\Box Sync\misclassifiedNetwork\data\';
filepath = 'C:\Users\xt9\Dropbox\Arthur_Xi_Xun\unobserved_links\final submission files\Data\';
load([filepath 'classDataV2.mat']);  clear majorProp;
global varIndiv varClass indY dobwt LB UB meanVI meanVC meanY options

%% SUMMARY OF VARIABLES
% Dependent variable: 
% - - math computation scale score (g3mathcomputss); or
% - - total math scale score (g3tmathss)
% Student demographics: 
% - - days absent (g3absent);
% - - motivation score (g3motivraw);
% - - grade 2 math raw score (g2mathsbsraw); 
% Class/teacher char.: 
% - - years of experience of teacher; 
% - - class type (g3classty)
varIndiv    = [10 11 12]; % individual chara
varClass    = 7; % 7; % class/teacher characteristic
indY        = 14;  % choice of outcome variables in data, either 13 or 14
dobCutoff   = 6.5; % 6.5
dobwt  = sum(dobDistrn < dobCutoff)/length(dobDistrn);
sizeCutoff  = 20;
options     = optimset('TolX',1e-8,'TolFun',1e-8,'MaxFunEvals',1e8,...
                       'MaxIter',1e8,'Display','off');

%% data management
for i = 1:2 % loop in the dummy 1('dobDistrn'<=6)
    for j = 1:2 % loop in small vs large classes
        %% partition the sample based on 'classSize' and 1('dobDistrn'<=6)
        index = find( (classSize > sizeCutoff ) == j-1 & (dobDistrn > dobCutoff ) == i-1 ); 
        % inpute characteristics of 'pseudo-members' in the partitioned
        % sample and demean the variables using conditional expection given group size
        [Yt{i,j},clsXt{i,j},indXt{i,j},LB(i,j),UB(i,j),...
            meanVI{i,j},meanVC{i,j},meanY{i,j},...
                OYt{i,j},OclsXt{i,j},OindXt{i,j}] = ...
             manage_STAR_data(classData(index),classSize(index)); 
    end;
end;

%% estimate individual, contextual, peer effects
inds = [1:15]; el = length(inds);
[BIND,BCLS,BCON,rBCON,rcm,rcmInd,rcmCls] = estSTAR(Yt,clsXt,indXt);
% save workspace for counterfactual analysis
save([filepath 'estSTAR.mat'],...
      'BIND','BCLS','BCON','rBCON','rcm','rcmInd','rcmCls',...
      'Yt','indXt','clsXt','meanVI','meanVC','meanY',...
      'OYt','OindXt','OclsXt','LB','UB');

%% bootstrap standard error 
B = 200; b = 1; count = 0;
BBIND = zeros(7,B); BBCLS = zeros(2,B); BBCON = zeros(2,B); 
BrBCON = zeros(2,B); Brcm = zeros(2,2,3,15,15,B);
while b < B+1, b
    % draw bootstrap sample
    for i=1:2, for j = 1:2,
         L = size(Yt{i,j},2);  bindex = randsample(L,L,1);             
         BYt{i,j} = Yt{i,j}(:,bindex);
         BindXt{i,j} = indXt{i,j}(:,:,bindex);
         BclsXt{i,j} = clsXt{i,j}(bindex);
    end; end;    
    [output1,output2,output3,output4,output5] = estSTAR(BYt,BclsXt,BindXt);
    if min(output1(1:2)) < .55 || max(output1(1:2)) > .99,
       count = count + 1; continue;
    else ...
        BBIND(:,b) = output1;
        BBCLS(:,b) = output2;
        BBCON(:,b) = output3;
        BrBCON(:,b)= output4;
        Brcm(:,:,:,:,:,b) = output5;
        b = b + 1;
    end;    
end;
save([filepath 'BSestSTAR.mat'],'BBIND','BBCLS','BBCON','BrBCON','Brcm');
% point est and confidence interval
[BIND, std(BBIND')'],

%% test: equal lambda for small & large classes
% (BIND(2)-BIND(1))/(std(BBIND(2,:)'-BBIND(1,:)')/sqrt(LB(1))), 
% [BIND(2)-BIND(1), std(BBIND(2,:)'-BBIND(1,:)')/sqrt(LB(1))], 

%% test: social interaction (all indivdidual linked with equal weights)
% under the null, row sum and column sum of reduced form coefficients
% should be the same
k = 3; inds = [1:6]; el = length(inds); 
for i = 1:2, 
    for j=1:2,                
        TM1 = diag(squeeze(rcm(i,j,k,inds,inds)/10)); 
        TM2 = reshape(squeeze(rcm(i,j,k,inds,inds)/10).*(ones(el,el)-eye(el)),el^2,1);
        TM2 = TM2(TM2~=0);
        RCM1 = []; RCM2 = [];
        for b = 1:B,
            tm = squeeze(Brcm(i,j,k,inds,inds,b)/10);
            RCM1 = [RCM1; diag(tm)'];
            tm2  = reshape(tm.*(ones(el,el)-eye(el)),el^2,1);
            RCM2 = [RCM2; tm2(tm2~=0)'];
        end;
        cov1 = cov(RCM1); cov2 = cov(RCM2); EL = length(TM2);
        R1 = [eye(el-1) zeros(el-1,1)] - [zeros(el-1,1) eye(el-1)];  
        R2 = [eye(EL-1) zeros(EL-1,1)] - [zeros(EL-1,1) eye(EL-1)];          
        wald1(i,j) =  (TM1(1:end-1)-TM1(2:end))'*inv(R1*cov1*R1')*(TM1(1:end-1)-TM1(2:end));
        wald2(i,j) =  (TM2(1:end-1)-TM2(2:end))'*inv(R2*cov2*R2')*(TM2(1:end-1)-TM2(2:end));
    end;
end;
chi2inv(.95,EL-1)
wald2,
        